Modelos Lineales

2019-02-11

Temario

Introducción
Regresión Lineal
- Simple
- Múltiple
Estudios observacionales
- Transversales
- Caso-Control
Modelos Lineales Generalizados
- Odds Ratio
- Prevalence Ratio
Estudios adicionales
- Cohortes - Risk Ratio
- Tiempo a evento - Hazard Ratio
- Desenlace conteo o tasa - Incidence Risk Ratio

Introducción

¿qué es R?

Software Libre

¿por qué usar R?

Ciencia Reproducible:
- limpiar bases de datos
- visualizar
- ejecutar modelos
- comunicar resultados

¿cómo usar R?

IDE: RStudio
Paquetes
Funciones

Introducción

¿por qué usar una regresión?

permite fijar al menos una variable independiente o exposición y observar una respuesta en la variable dependiente o desenlace.
permite explicar el cambio promedio de un evento \(Y\) en base a cambios en \(X\), usando coeficientes o medidas de asociación.
permite predecir la probabilidad asociada a un evento.
.

Regresión Lineal Simple

características

una variable independiente (simple)
una variable dependiente (univariada)
ambas variables deben ser numéricas

objetivos

ajustar datos a una recta
interpretar medida de bondad de ajuste ( R² ) y coeficientes
evaluar supuestos
visualizar el modelo

\[ Y = \beta_0 + \beta_1 X_1 + \epsilon \]

RLS: datos

espir <- read_dta("data-raw/espirometria.dta") %>% as_factor()
espir

## # A tibble: 654 x 7
##     caso codigo  edad   vef talla sexo   fumar
##    <dbl>  <dbl> <dbl> <dbl> <dbl> <fct>  <fct>
##  1     1    301     9  1.71  57   female No   
##  2     2    451     8  1.72  67.5 female No   
##  3     3    501     7  1.72  54.5 female No   
##  4     4    642     9  1.56  53   male   No   
##  5     5    901     9  1.89  57   male   No   
##  6     6   1701     8  2.34  61   female No   
##  7     7   1752     6  1.92  58   female No   
##  8     8   1753     6  1.41  56   female No   
##  9     9   1901     8  1.99  58.5 female No   
## 10    10   1951     9  1.94  60   female No   
## # ... with 644 more rows

RLS: distribución

supuestos:
- #1 linealidad
- #2 independencia de observaciones

espir %>% 
  ggplot(aes(edad,vef)) +
  geom_point()

error = residual

RLS: suma de mínimos cuadrados

Cálculo de la sumatoria del cuadrado de los residuales hacia la media y la recta:

\[ SSE(mean) = \sum (data - mean)^2 \] \[ SSE(fit) = \sum (data - fit)^2 \]

\[ Var(x) = \frac{SSE(x)^2}{n} \]

Medida de bondad de ajuste:

\[ R^2 = \frac{Var(mean) - Var(fit)}{Var(mean)} \]

RLS: R²

wm1 <- lm(vef ~ edad, data = espir)
wm1 %>% glance() %>% select(r.squared:df)

## # A tibble: 1 x 6
##   r.squared adj.r.squared sigma statistic   p.value    df
## *     <dbl>         <dbl> <dbl>     <dbl>     <dbl> <int>
## 1     0.572         0.572 0.568      872. 2.45e-122     2

INTERPRETACIÓN

edad explica el 57% de la variabilidad de VEF
existe un 57% de reducción en la variabilidad de VEF al tomar en cuenta la edad

RLS: coeficientes

wm1 %>% tidy()

## # A tibble: 2 x 5
##   term        estimate std.error statistic   p.value
##   <chr>          <dbl>     <dbl>     <dbl>     <dbl>
## 1 (Intercept)    0.432   0.0779       5.54 4.36e-  8
## 2 edad           0.222   0.00752     29.5  2.45e-122

wm1 %>% confint_tidy()

## # A tibble: 2 x 2
##   conf.low conf.high
##      <dbl>     <dbl>
## 1    0.279     0.585
## 2    0.207     0.237

\(\beta_{edad}\): En la población,
por cada incremento de edad en una unidad, el VEF en promedio incrementa en 0.22 mL/s,
con un intervalo de confianza AL 95% de 0.21 a 0.24 mL/s.
Este resultado es estadísticamente significativo con un valor p < 0.001

RLS: supuesto #3 normalidad

wm1 %>% 
  augment() %>% 
  ggplot() +
  geom_qq(aes(sample=.std.resid)) + 
  geom_qq_line(aes(sample=.std.resid))

META: todos los puntos sobre la línea

wm1 %>% augment() %>% 
  ggplot(aes(.std.resid)) +
  geom_histogram()

RLS: supuesto #4 homoscedasticidad

META: distribución idéntica a ambos lados de la línea

wm1 %>% 
  augment() %>% 
  ggplot(aes(.fitted,.std.resid)) +
  geom_point() +
  geom_smooth() +
  geom_hline(yintercept = c(0))

RLS: ¿cómo se ve el modelo?

\[ Y = \beta_0 + \beta_1 X_1 + \epsilon \] \[ VEF = 0.60 + 0.22 (edad) + \epsilon \]

espir %>% 
  ggplot(aes(edad,vef)) +
  geom_point() + 
  geom_smooth(method = "lm")

RLS: retroalimentación

R² indica el porcentaje de variabilidad del desenlace (var. dependiente) explicada por el predictor (var. independiente).
los coeficientes permiten predecir el posible valor del desenlace en base a un modelo estadístico.
los supuestos permiten evaluar qué tan adecuado es el ajuste de los datos al modelo.

RLS: ejercicio

ajustar una recta

espir %>% 
  ggplot(aes(edad,talla)) +
  geom____()

identificar coeficientes y R²

# recordar: y ~ x
wm1 <- lm(_____ ~ _____, data = _____)
wm1 %>% g_____
wm1 %>% t_____
wm1 %>% c_____

evaluar supuestos: normalidad y homoscedasticidad

wm1 %>% augment() %>% 
  ggplot() + _____

wm1 %>% augment() %>% 
  ggplot(aes(.fitted,.std.resid)) + __________

Intermedio #1: Tabla 1 y 2

“… ¡en un solo comando!”

compareGroups(fumar ~ edad + vef + talla + sexo,
              data = espir, byrow=T#,method=c(vef=2)
              ) %>% 
  createTable(show.all = T) #%>% export2xls("table/tab1.xls")

## 
## --------Summary descriptives table by 'Habitos de fumar'---------
## 
## ________________________________________________________ 
##               [ALL]        No          Si      p.overall 
##               N=654       N=589       N=65               
## ¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯ 
## edad       9.93 (2.95) 9.53 (2.74) 13.5 (2.34)  <0.001   
## vef        2.64 (0.87) 2.57 (0.85) 3.28 (0.75)  <0.001   
## talla      61.1 (5.70) 60.6 (5.67) 66.0 (3.19)  <0.001   
## sexo:                                            0.071   
##     male   336 (51.4%) 310 (92.3%) 26 (7.74%)            
##     female 318 (48.6%) 279 (87.7%) 39 (12.3%)            
## ¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯

Intermedio: distribuciones

espir %>% 
  select(edad,vef,talla) %>% 
  gather(key,value) %>% 
  ggplot(aes(value)) +
  geom_histogram() +
  facet_wrap(~key,scales = "free")

## # A tibble: 1 x 10
##   n_obs   min   max  mean    sd   q25   q50   q75 skewness kurtosis
##   <int> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>    <dbl>    <dbl>
## 1   654 0.791  5.79  2.64 0.867  1.98  2.55  3.12    0.661     3.30

Regresión Lineal Múltiple

características

dos o más variables independientes (múltiple)
una variable dependiente numérica

objetivo

epidemiológico: controlar confusión en el análisis.
- estandarización (directa/indirecta),
- estratificación,
- puntajes de propensión,
- regresión multivariable
estadístico: emplear múltiples predictores
- lineales, no-lineales, segmentadas

RLM: distribución

espir %>% 
  ggplot(aes(edad,vef,colour=sexo)) +
  geom_point()

RLM: R² y coeficientes

## # A tibble: 1 x 6
##   r.squared adj.r.squared sigma statistic   p.value    df
## *     <dbl>         <dbl> <dbl>     <dbl>     <dbl> <int>
## 1     0.607         0.606 0.544      503. 9.49e-133     3

edad y sexo explica el 60% de la variabilidad de VEF.

## # A tibble: 3 x 7
##   term        estimate std.error statistic   p.value conf.low conf.high
##   <chr>          <dbl>     <dbl>     <dbl>     <dbl>    <dbl>     <dbl>
## 1 (Intercept)    0.605   0.0781       7.74 3.79e- 14    0.451     0.758
## 2 edad           0.220   0.00722     30.6  7.32e-128    0.206     0.235
## 3 sexofemale    -0.323   0.0426      -7.59 1.13e- 13   -0.407    -0.240

\(\beta_{sex:female}\)
En la población, el VEF en promedio en mujeres es 0.32 mL/s menor que en hombres,
con un intervalo de confianza AL 95% de 0.24 a 0.41 mL/s, ajustando por edad.
Este resultado es estadísticamente significativo con un valor p < 0.001

RLM: supuesto normalidad y homoscedasticidad

wm1 %>% 
  augment() %>% 
  ggplot() +
  geom_qq(aes(sample=.std.resid)) + 
  geom_qq_line(aes(sample=.std.resid))

wm1 %>% 
  augment() %>% 
  ggplot(aes(.fitted,.std.resid)) +
  geom_point() +
  geom_smooth() +
  geom_hline(yintercept = c(0))

RLM: post-estimación

multicolinealidad: VIF
outliers: puntos influenciales o ruido

RLM: ¿cómo se ve el modelo?

\[ Y = \beta_0 + \beta_1 X_1 + \beta_2 X_2 + \epsilon \] \[ VEF = 0.60 + 0.22 (edad) - 0.32 (sexo) + \epsilon \]

wm1 %>% augment() %>% 
  ggplot(aes(colour= sexo)) +
  geom_point(aes(edad,vef),alpha=0.05) +
  geom_line(aes(edad,.fitted), lwd=1.5)

RLM: ejercicio

3 predictores

wm1 <- lm(vef ~ edad + sexo + fumar, data = espir)

R² y coeficientes

wm1 %>% g___()
wm1 %>% t___() 
wm1 %>% c___()

supuesto normalidad y homoscedasticidad

wm1 %>% augment() %>% 
  ggplot() +
  geom_qq________ + 
  geom________

wm1 %>% augment() %>% 
  ggplot(aes(_____,_____)) +
  geom_p____()

Intermedio #2: Modificación de efecto

¿cómo se ve el modelo?

\[ Y = \beta_0 + \beta_1 X_1 + \beta_2 X_2 + \beta_3 (X_1 * X_2) + \epsilon \]

## # A tibble: 4 x 4
##   term         estimate conf.low conf.high
##   <chr>           <dbl>    <dbl>     <dbl>
## 1 (Intercept)     0.253   0.0911     0.416
## 2 edad            0.243   0.226      0.259
## 3 fumarSi         1.94    1.13       2.76 
## 4 edad:fumarSi   -0.163  -0.223     -0.102

\(\beta_{interacción}\)
La diferencia en la media de VEF asociado al incremento en 1 año de edad en fumadores
menos
la diferencia en la media de VEF asociado al incremento en 1 año de edad en no fumadores
es -0.16 mL/s, con un IC 95% de -0.22 a -0.1.

Estudios Observacionales

Asociar:
- condición (expuesto o no-expuesto) y
- resultado (caso o control)
El tipo de análisis depende del tipo de estudio:
- “Reclutamos \(x_1\) sujetos expuestos y \(x_2\) sujetos no expuestos”
- Tipo Cohorte
- “Reclutamos \(x_1\) sujetos caso y \(x_2\) sujetos control”
- Tipo Caso-Control

tipo de estudio	medida de asociación
caso - control	OR
cohorte	RR
transversal	PR
tiempo a evento	HR

Bases de datos

Cohorte

riskdb <- read_dta("data-raw/epidb.dta") %>% as_factor() %>% 
  mutate(aids=if_else(aids==0,"No","Yes") %>% as.factor(),
         ccr2=if_else(ccr2==1,"No","Yes") %>% as.factor())

Encuesta a nivel nacional

logdb <- read_dta("data-raw/brfss2000.dta") %>% as_factor() %>% 
  filter(!is.na(depress))
  #mutate(dn_cro=as.factor(dn_cro))

Cohorte

Elegimos los totales de expuestos y no expuestos

riskdb %>% 
  count(ccr2,aids) %>% 
  spread(aids,n) %>% 
  mutate(risk=Yes/(No+Yes)) %>% 
  mutate(rr=.$risk[2]/.$risk[1])

## # A tibble: 2 x 5
##   ccr2     No   Yes  risk    rr
##   <fct> <int> <int> <dbl> <dbl>
## 1 No      222   149 0.402 0.821
## 2 Yes      63    31 0.330 0.821

Si cambiamos el tamaño del grupo (no-expuestos x 100), no cambia el Riesgo Relativo (RR)

## # A tibble: 2 x 5
##   ccr2     No   Yes  risk    rr
##   <fct> <dbl> <dbl> <dbl> <dbl>
## 1 No    22200 14900 0.402 0.821
## 2 Yes      63    31 0.330 0.821

Caso Control

Elegimos los totales de casos y controles
Si cambiamos el tamaño del grupo (casos x 100), sí cambia el Riesgo Relativo (RR)

riskdb %>% 
  count(ccr2,aids) %>% 
  spread(aids,n) %>% 
  mutate(Yes=Yes*100,
         risk=Yes/(No+Yes)) %>% 
  mutate(rr=.$risk[2]/.$risk[1])

En estudios caso-control usamos el Odds Ratio (OR)

## # A tibble: 2 x 7
##   ccr2     No   Yes  risk  odds    rr    or
##   <fct> <int> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 No      222 14900 0.985  67.1 0.995 0.733
## 2 Yes      63  3100 0.980  49.2 0.995 0.733

sobreestima el Riesgo Relativo (RR)